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Abstract 

The Hubbard model is studied in which disorder is introduced by putting 
the on-site interaction to zero on a fraction / of (impurity) sites of a square 
lattice. Using Quantum Monte Carlo methods and Dynamical Mean-Field 
Theory we find that antiferromagnetic long-range order is initially enhanced 
at half-filling and stabilized off half- filling by the disorder. The Mott-Hubbard 
charge gap of the pure system is broken up in two pieces by the disorder: one 
incompressible state remains at average density n = 1 and another can be 
seen slightly below n = 1 + /. Qualitative explanations are provided. 
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1 Introduction 



The problem of the interplay between disorder and interactions in systems of elec- 
trons is challenging and has a long history . Disorder can, by itself, localize elec- 
trons and make the conductivity disappear, the Anderson transition. At appropriate 
densities, interactions can also cause the formation of insulating states via the Mott 
transition. In the latter case, the insulating states often possess additional magnetic 
or charge ordering. The simultaneous presence of disorder and interactions is re- 
alized in many systems in nature, with fundamental and fascinating consequences, 
e.g. the unusual magnetic properties of heavily doped, compensated semiconductors, 
and the unique transport properties of two-dimensional and bulk superconductors 
[H]. If both disorder and interactions are present, the effects can sometimes reinforce 
each other, but can also compete. In classical spin models, it has recently become 
clear that quenched disorder can turn (temperature-driven) symmetry-breaking first 
order phase transitions into continuous phase transitions, possibly with new inter- 
vening critical points 0. Although the arguments do not appear to apply directly 
to interacting fermions, disorder may lead to similarly interesting phenomena. 

The Hubbard model for interacting electrons on a lattice 0| exhibits both Mott 
metal-insulator and magnetic phase transitions. On the one hand, it is expected 
that the interaction U induces a gap at half-filling by separating many-body states 
with doubly occupied orbitals from those with holes, when the on-site repulsion 
becomes larger than the non-interacting bandwidth. On the other hand, also near 
half-filling, there is a tendency to antiferromagnetic (AF) ordering of the electron 
spins. One could now ask the question whether the Mott transition in fact ever 
occurs in the absence of some associated symmetry breaking such as magnetic order 
[^]. For interacting bosons the Mott transition (in this case from a superfiuid) to an 
insulator is known to occur without any other accompanying order parameter 0. 

To investigate the behavior of disordered interacting electron systems it is of 
interest to consider the Hubbard Hamiltonian into which disorder is introduced: 

where Cjo- is the annihilation operator for an electron at site j with spin a. tij 
is the hopping integral (non-zero only between neighboring sites i and j), Uj the 
on-site interaction at site j between electrons of opposite spin, and n^o- = cJ^Cj^. 
is the occupation number operator. Various kinds of disorder are allowed for in 
(IID: (i) in the hopping integrals tij, (ii) in the on-site energy ej, or (iii) in the on- 
site interaction Uj. Traditionally, because the case of non- interacting electrons was 
considered first, attention has been devoted to the first two kinds of disorder, termed 
off- diagonal and diagonal disorder, respectively. In recent Quantum Monte Carlo 
simulations, the first two kinds were studied including uniform interactions, showing 
the destruction of AF long-range order at half-filling for some critical strength of 
disorder in case (i), and allowing for few conclusions in case (ii) because of sign 
problems already at half-filling p. In the present paper, we study case (iii), i.e. the 
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Hamiltonian: 



H=-t clcj, + Uj {rij^ - \){nj^ - i) - ^ ^ n^,, (2) 

in which the on-site repulsion Uj is turned off from on a fraction / of all sites ((i, j) 
denotes neighboring sites i and j), i.e. the Uj are taken from a bimodal distribution, 

P{U,) = {l-f)5{U,~U) + f5{U,) . (3) 

A physical realization of this kind of disorder could be impurity atoms in a crystal 
which can accomodate two electrons with a strongly reduced Coulomb repulsion, 
e.g. non-magnetic impurities (like Zn) replacing copper atoms in the copper oxides. 
Here, we do not want to attach ourselves to a specific physical system, but try to 
understand what happens to ordered states of interacting electrons when disorder 
of kind (iii) is introduced. One expects the Mott transition to be shifted to average 
density n = 1 + / (at which density also the U ^ Q sites start to be doubly occupied), 
so that it is perhaps separated from the AF magnetic order which is likely to stay 
at n = 1. New types of ordered states may also result. 

To explore these phenomena tractable and reliable theoretical methods are hard 
to come by. Specifically, mean-field approaches to the physics of strongly corre- 
lated quantum systems, like the Hartree-Fock approximation and the slave-boson 
mean-field approach, can provide significant insight into the possible phases which 
can arise due to interaction effects, but are well-known to overestimate significantly 
the tendency to form ordered states [0, |^. Other analytical approaches like the 
renormalization group, while powerful, likewise have their limitations [0]. The most 
obvious numerical approach, exact diagonalization of the Hamiltonian, is restricted 
to small lattices (10-20 particles) 0. In this paper, we intend to study the ef- 
fects of interactions and randomness by means of a stochastic sampling of the state 
space, the Quantum Monte Carlo (QMC) method. In this approach, correlation 
effects are treated exactly. A few hundred interacting electrons can be simulated, 
roughly an order of magnitude larger than by exact diagonalization. While the 
constraints of this finite size are still very significant, this number of particles is 
often sufficient to allow systematic extrapolations to bulk behavior, at least in two 



dimensions. The QMC method we use is the Determinant Monte Carlo method flO 
which has been extensively applied to non-random interacting fermion Hamiltoni- 
ans. We combine these simulations with calculations using Dynamical Mean-Field 



Theory (DMFT), or: the infinite-dimension approach. |TT|. This technique improves 
significantly upon the Hartree-Fock approximation since it is non-perturbative and 
takes quantum fluctuations into account via a dynamical self-energy. The infinite- 
dimension approach allows us to analyze problems in the thermodynamic limit and 
in larger ranges of parameter space, e.g. lower temperatures, than accessible to 
the determinant algorithm. DMFT has provided valuable insight into the physics 
of strongly correlated electron systems, e.g. the Mott-Hubbard transition 0, [Tl[| . 
Both techniques, which recently have been applied to disordered Hubbard models 
also P, |12|, will be briefly described in Section 2. In subsequent sections, we 
introduce the quantities of interest that are calculated, present results and discuss 
them. 
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2 Quantum Monte Carlo and Dynamical Mean 
Field Theory 

2.1 Determinant Monte Carlo method 

A brief description of the Determinant Monte Carlo algorithm for the Hubbard 
Hamiltonian (|^) is given. One is interested in computing operator expectation values 
like: 

Tr Ae'f^^ 

- ■ w 

The trace over the fermion degrees of freedom cannot be performed analytically 
due to the quartic interaction term. In order to reduce the problem to a quadratic 
Hamiltonian we discretize, (3 = L^-At, and employ the Trotter decomposition |T^. 
The partition function can now be expressed as: 

Z = Ti e-""" = Tr [e"^"^]^- ^ Tr [e-^"^e'^"^]^- . (5) 

Here K includes the a priori quadratic pieces of H, the hopping and chemical po- 
tential terms, while P is the on-site interaction term. In the second step, errors in 
measured quantities of order t?7(Ar)^ are introduced. Now that P has been iso- 



lated, the discrete Hubbard- Stratonovich transformation of Hirsch [T^ can be used 
to decouple the interaction. If U > 0, 

Here cosh(AAr) = exp([/Ar/2). A field variable Si^r (of Ising type) must be intro- 
duced at each lattice site and imaginary time slice. If ?7 < the field Si^r couples to 
the charge rii^ + riii on each site. In either case, the end result is that the partition 
function can be rewritten as: 

Z = 5]Trne-V^'(') , (7) 

where the fermion operators appear only quadratically in the transformed P' and 
/ denotes the Trotter time slice. The interacting fermions have been replaced by 
non-interacting fermions evolving in a space- and time-dependent field. Performing 
the fermion trace is now possible, and results in: 

Z = J2 detM^{Si,r) detMi{S^,r) , (8) 

where the detailed forms of the matrices M„, which have dimension the number of 
spatial lattice sites N, are written down by Blankenbecler et a/.|]lO|. 

One crucial bottleneck in the algorithm is that while at high temperatures the 
product of the fermion determinants is positive and can be used as a statistical 
weight, at low temperatures (high (3) when U > and the system is away from 
half-filling, this is no longer always true. In this case, as (3 increases the product 
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can become negative. This is referred to as the "sign problem" . The sign problem 
precludes simulations for large lattices at low temperatures. When [/ < no sign 
problem occurs, and this is one reason the attractive Hubbard model is easier to 
study. 

The Monte Carlo simulation proceeds by going through the space-time lattice of 
Si^r and computing the ratio of the product of fermion determinants before and after 
a single spin is flipped. Naively, the simulation time scales as the fourth power of 
the number of sites A^. One power of comes from updating all the sites. For each 
such update, the calculation of the determinant is of order A^^. However, one can 
reformulate the algorithm to scale with one lower power of as follows: instead of 
evaluating the new determinant from scratch, the ratio of the determinants may be 
evaluated in terms of the elements of M~^, the fermion Green function. One then 
must keep track of M~^, which of course changes as moves are accepted. However, 
the update of resulting from a change to Mg. is only an A^^ procedure, due to 
the specific local nature of the change in M„. In summary, the technique involves 
the manipulation (mostly multiplication) of matrices of dimension the number of 
spatial sites in the system. This number may be up to several hundred on present- 
day computers. 

The inputs to the computer code are the quantities that specify the physical 
system. They are: the spatial size of the system and the inverse temperature (3. Also, 
the parameters in Hamiltonian (Q), t,Uj, and n, need to be specified. Finally, the 
parameters of the simulation, the Trotter time interval At (which then determines 
the number of time steps Lr), the number of warm-up and measurement sweeps, 
the time between measurements etc. must be selected. Quite some care needs to 
be taken in determining these latter quantities to ensure the simulation results are 
meaningful. 

One crucial feature of the Determinant Monte Carlo simulations is that the 
fermion Green functions Gij{T,T') = {Tci{r)c^j{r')) are readily obtainable. In fact, 
G is just the inverse of the matrix M whose determinant gives the Boltzmann 
weight. It is available "free of charge" (at least for r = r') in the sense that it 
has already been constructed in moving the field variables. With the single-particle 
Green function in hand, we can measure quantities involving more than two fermion 
operators by performing the appropriate Wick contractions and expressing them in 
terms of products of single-particle propagators. There is considerable flexibility in 
the quantities one can calculate, and it is possible to decompose diagrams and sep- 
arate self-energy and vertex contributions. This is useful when comparing quantum 



simulations with analytic calculations 16 



The Determinant Monte Carlo algorithm as described above has been extensively 
used for over a decade now. The more reliable applications have mainly been due to 
a number of recent technical advances. Here we only mention the most important 
one for our present studies: the introduction of matrix stabilization procedures has 



significantly improved the ability to reach low temperatures |]17| , p!8| . With these 
techniques inverse temperatures as high as f3t = 20 can be attained. 

In the present study, we have used U = 8t (equal to the band width) combined 
with a Trotter time step At = l/12t, ensuring small Trotter errors. Experience 



from previous studies |18] and independent tests in the present work teach that 
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an inverse temperature /3t = 8 is sufficiently high to effectively yield ground-state 
thermodynamic and correlation functions on the lattices considered. This tempera- 
ture is used throughout. We typically use 200 to 300 warm-up sweeps and 2000 to 
4000 measurement sweeps. A simplifying detail of the present problem, in which a 
fraction / of the sites have U = 0, is that the corresponding Hubbard-Stratonovich 
coupling A in @ equals zero and the local electron Green function does not need 
updating. We study square lattices with linear size up to L = 10 and average over 
up to 40 disorder realizations. 

Like the original Hubbard model, Hamiltonian preserves a particle-hole sym- 
metry at half-filling [n = 1 and fi = 0), i.e. it is invariant under the "staggered" 
particle-hole transformation cj^ — > (— l)*Cjg.. This symmetry ensures that there is 
no sign problem at n = 1. Note that the particle-hole symmetry corresponds to 
different chemical potentials on the two constituents. When going off half-filling we 
do encounter the sign problem. However, we find it to be less detrimental than in 
the case without disorder. Aided furthermore by the disorder averaging it is possible 
to obtain meaningful data off half-ffiling (see also section ^). We further support 
our findings by comparison to results of DMFT calculations. 



2.2 Dynamical Mean— Field Theory 

For classical spin models (e.g. the Ising model) it is well known that the Weiss 
molecular field (or: mean field) theory becomes exact in the limit of high spatial 
dimensions. For lattice electrons this limit was introduced only recently; with the 
proper scaling of the hopping element in (|^), t = t*/v^ {z is the number of nearest 
neighbors) it leads to a quantum mechanical dynamical mean-field theory (DMFT). 
Here we give a brief description of this approach to our problem at hand; more 



details on the method are found in the literature |[TT|, |T3 . 

Thermodynamic properties are determined from the averaged grand potential 

/Jfiav = -(( In Tr exp{-0))),, , (9) 

where (. . .)av denotes the disorder average. The main simplification in the limit 
d ^ oo lies in the reduction of the average, which has to be performed only on a 
single site |]T^ (for each allowed value of Uj, and using a coherent potential type of 



approximation). The explicit expression for the grand potential in the paramagnetic 
phase is given by: 

/oo 

+ Nj2\nG-^-N{\nZ{G,^,U,}),^, (10) 

cm 

where uJn = {2n+l)7T / j3 are Matsubara frequencies and N^{E) is the density of states 
(DOS) of non-interacting electrons (scaled to be suitable for comparison with the 
square lattice [^). In (|10D also appear: the generalized atomic partition function Z 



(shown to be equivalent to that of a single-impurity Anderson model p2|]), and the 
complex quantities Gan and which at this point enter as variational parameters. 
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The physical values of Gan and So-n (namely, the local Green function, Gan = Gu^an, 
and the electron self-energy. So-™ = Sjj respectively) correspond to those for which 
( p!OD is stationary |21|]. The stationarity conditions: 6flav/SGun = , 6Qav/S^an = 
0, yield two coupled sets of self-consistent equations for G^n and 



-'an- 



oo 



G.n = - I dE-—^^ (11) 



-oo 

/3 



Gan = J dr e^--^{cAT)c:mT).. (12) 



where {■ ■ .)t denotes the thermal average (which depends on G„n and So-n implic- 
itly). Physically, the Dyson equation (|11|) describes the local Green function of 
independent electrons moving in a homogeneous dynamical potential So-™. It can 
be solved by a simple integration for each Matsubara frequency. The functional 
integral (|l^) on the other hand is highly non-trivial since it couples all Matsubara 
frequencies. The interacting, disordered problem is mapped onto an ensemble of 
single-impurity problems, complemented by a self-consistency condition which in- 
troduces the lattice into the problem |T2| , |19[. We emphasize that, in contrast to 
conventional mean-field theories, in this dynamical mean-field theory (DMFT), the 
action remains time dependent, i.e. local fluctuations are retained. 

The local interacting problem ([T^ ) is solved numerically using a finite-temperature, 
auxiliary-field Quantum Monte Carlo (QMC) method |2^. Quite similar as in the 



finite-dimensional QMC algorithm described above, the electron-electron interac- 
tion is formally replaced by an interaction of independent electrons with a dynamical, 
auxiliary field of Ising-type spins. To this end the interval [0,/3] is again discretized 
into Lt- steps of size At = P/Lt-. Equivalently, there is a high energy cut-off of 
Matsubara frequencies, i.e. = \2n + l|7r//5 < tt/At, n = —Lt-/2, . . . ,Lt-/2 — 1. 
The computer time grows like oc f3^, restricting to values below ~ 150 and 
(3t* < 50 ... 70 on present supercomputers. For small Lr {Lr < 20) one can perform 
a full enumeration (instead of Monte Carlo sampling) of all 2^^ possible configu- 
rations of the auxiliary field. Sign problems turn out to be absent in the DMFT 
approach. 

Self-consistency is obtained iteratively as follows: the Green function G (omit- 
ting indices) is calculated from some initial self energy, e.g. S = 0, by the Dyson 
equation ([TT|) . Then, the new Green function Gnew is determined by solving (|12D for 
each value of Uj with the QMC method, and averaging. Finally, the calculation of 
the new self energy Snew = S — G^^^ + G^^ completes one iteration. Typically 10 
iterations with 20000 MC sweeps are necessary to obtain a convergence of ~ 10^'^. 
Close to a phase transition the convergence is much slower and the statistical errors 
are larger due to strong fluctuations. At large [/-values {U > 12t) the Monte Carlo 
sampling becomes more and more inefficient due to "sticking" problems, i.e. there 
are two (or more) minima in the free energy and the single spin-flip algorithm is no 
longer able to transfer between them. 

The self-consistent equations ([TT|) and (|5) can easily be extended to a phase 



with AF long range order (For details of the implementation see [0). One uses the 
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property that, in the case of finite sublattice magnetization but zero total magnetiza- 
tion, i.e. excluding ferrimagnetism, the Green function for an f —spin on sublattice 
A is equal to that of a | —spin on sublattice B and vice versa. 



3 Quantities of Interest 

In this paper, we only study thermodynamic quantities and static correlation func- 
tions. At a later stage also (imaginary-) time-dependent correlation functions will 
be of interest, e.g. to study frequency-dependent response functions. 

In the QMC calculations, the (equal-time) electron Green function Gij^ = 
(c^^c^^) = Mj~^, described in section immediately leads to electron densities. 



(^icr) = 1 ~ Giia, sjid statlc Correlation functions. The magnetic correlations which 
give crucial information on the physics of the Hubbard model can be monitored by 
measuring the equal-time spin-spin correlation function c(/) = {{uji — nj|)(nj+/| — 
rij+il)) and its Fourier transform S{q) = I]; e*'' 'c(/), the magnetic structure factor. 
Specifically, a finite-size scahng study of ^(Tr, vr) leads to the sublattice magneti- 
zation M of the phase with AF long range order, according to spin wave theory 
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%^^fH-0(l/L), (13) 

where L is the linear system size. 

In the DMFT calculations (see section p.2|) , the resulting Green functions also 



lead to, for instance, the densities no- and sublattice magnetization M. The Green 
function Ga^j is just the Fourier transform of G^n in section p.2| . Note that in this ap- 
proach quantities like the order parameter M are calculated in the thermodynamic 
limit, so no finite-size scaling is required. The vanishing of M is used to determine 
phase boundaries as a function of parameters like /i, U, T, or /, leading to e.g. the 
Neel temperature Tn. A different way to obtain these phase boundaries is to cal- 
culate the corresponding susceptibility xaf (via two-particle correlation functions; 
see [O for details). Divergence of xaf signals a continuous phase transition. 



4 Results and Discussion 
4.1 Magnetic order 

First, we study the effect of an increasing concentration f oi U = sites on the 
stability of AF long range order at n = 1. The static AF structure factor 5(7?, vr) 
is calculated a.t U = 8t for different lattice sizes and temperatures. For L < 10, 
S^TT, vr) is found to saturate at T t/8. /^From the saturated values the ground-state 
sublattice magnetization M can be extrapolated using finite-size scaling, according 
to (0). Scaling plots for different f at U = 8t are shown in Fig. 1. For / < 0.36, 
S {it , it) / L"^ extrapolates to a finite value in the thermodynamic limit. Long range 
order has disappeared for / > 0.5. Fig. 2 presents the extrapolated values M(/). For 
small disorder M is found to increase with /, it reaches a maximum around / = 0.1, 
and eventually vanishes at fc ~ 0.4. This initial increase of the order parameter 
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T=t/8 U = 8t 




Figure 1: Finite-size scaling of the antiferromagnetic structure factor ^(Tr, tt) on 
the square lattice at half-filling for / = 0, 0.125, 0.188, 0.25, 0.36, 0.5 (top to 
bottom; / is the fraction of f/ = sites), computed using the Determinant QMC 
method with U = 8t and T = t/8. L denotes the linear system size and N the 
number of sites (A^ = L^). For values of / that do not correspond to an integer 
number of defects, we have interpolated between the results for the two bracketing 
concentrations. 



prior to the formation of a disordered phase is remarkable. Similar phenomena have 
however been seen in studies of the Hubbard model with random site energies in 
infinite dimensions (an initial increase with disorder of the critical temperature for 
AF order [0). Local enhancement of AF fluctuations by vacancies has also been 



reported recently p5 



For the present model, at half-filling, we do not find initially enhanced AF order 
in the infinite-dimension approach as can be seen in Fig. 3 (Due to the neglect 
of spatial fluctuations, we do of course find that AF order survives longer when 
increasing disorder: fc ~ 0.75. This can be higher than the percolation threshold 
because no sites are actually removed). On the other hand, also in this approach 
we can see that AF order is enhanced by introducing U = sites: we find that AF 
order occurs at densities for which the pure system is not ordered. A calculated (/, n) 
phase diagram shows that AF order persists up to increasing values of n (relative to 
/ = 0, for which in DMFT AF order occurs up to n = 1.14, i.e. somewhat off half- 
filling) when the disorder increases from / = to about / = 0.36 |2^. Therefore, 



again the presence of impurity sites stabilizes the magnetic order. We note that 
the calculation of this phase diagram using the Determinant QMC method in finite 
dimension is very hard, if not impossible, because of the sign problem. We, however, 
expect the phase boundary between AF and disordered phase to be of similar form 
(but shifted towards half-filling). 
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U = 8t 




Figure 2: Ground-state staggered magnetization M as a function of fraction / of 
U — sites, as extrapolated from Fig. 1. 



U=8t 




Figure 3: Staggered magnetization M and Neel temperature as a function of 
fraction f oiU — sites, computed using Dynamical Mean Field Theory (DMFT). 
Tn values are obtained by extrapolating At to zero (see section 2.2) those for M 
are not, since extrapolation is cumbersome. Also, small values of M are hard to 
compute. 
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To come to an explanation for this disorder-stabilized magnetic order we argue 
that the situation is different in one important aspect from the case of electron- 
or hole-doping of the antiferromagnet at half-filling in the pure Hubbard model, 
which destroys AF order rapidly. In the pure case, the extra electrons (or: holes) 
are mobile and therefore particularly effective in destroying long-range order. With 
U = sites, however, extra particles will be localized on these sites, where they 
don't have to pay the on-site repulsion energy. So whereas magnetic moments are 
destroyed locally (on the U = sites), the long-range order originating from the 
U ^ sites will be able to persist much longer. That U = sites are not so 
detrimental to AF order can also be seen from the fact that the usual perturbative 
argument to obtain an AF (super exchange) Heisenberg coupling in the pure model 
still goes through with an intermediate U = site. Finally, a qualitative way to 
understand the enhanced M at n = 1 comes from spin wave theory: normally the 
quantum fluctuations of spin waves reduce the magnetization. The effect of disorder 
can be thought of as to introduce a lifetime, which can be seen to reduce the effect 
of the spin waves on M, thereby enhancing it. Analogously, introducing U = sites 
into the negative— f/ Hubbard model we expect to lead to enhanced charge-density 
wave order and superconductivity. The spin-wave argument may also explain the 
absence of the enhancement of M in DMFT in Fig. 3: in this approach the coherent- 
potential approximation (CPA) is used to treat disorder (see section [2.2[ ) and the 
CPA is known not to describe localization (in this case of spin waves), because it 
does not introduce spatial fluctuations iT3[| . 

4.2 Charge order 

A second type of order that is of interest is charge order, which generally is signalled 
by a range of constant density (or: zero compressibility k = dn/dii) when the 
chemical potential is varied. For instance, in the Quantum Hall effects the occurrence 
of incompressible states is the central theme. In the pure Hubbard model a charge 
gap (proportional to U) occurs around n = 1: the system resists against having 
doubly occupied sites. To study the effect of disorder by having U = sites, we 
calculate the (disorder-averaged) density n when varying /i in (|^) at f/ = 8t for a few 
values of impurity concentration /. To get a more detailed picture, we discriminate 
between U = and U = 8t sites. The results are displayed in Figs. 4 to 7. 

Before embarking on a discussion of the results, we first discuss the specific 
averaging procedure that has been used to obtain the result for n(/i) in Figs. 4 and 
5 (using Determinant QMC on 6 x 6 lattices). Because of the sign problem it is not 
possible to obtain reliable information at low enough temperatures on large lattices 
in the pure model for densities 1 < < 1.3 ([[l^], see also section |2.1|) . Introducing 
disorder turns out to help us in two ways: (i) in cases where the sign problem is less 
severe, e.g. for n = lA, the computed average sign is closer to 1 in the presence 
of disorder, which means that the problem is more tractable, and (ii) the disorder 
average that one has to perform leads to a natural way of taming the sign problem. 
Concerning (ii): since for each realization of (quenched) disorder a relatively small 
part of phase space (but hopefully large enough to be representative) is sampled, it 
can happen that for some realizations one is hit by the sign problem, but for others 
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T=t/8 U=8t 6x6 




1 2 



Figure 4: Average total density n as a function of jj, for different impurity concen- 
trations / = 0, 0.056, 0.111. Dotted lines indicate the corresponding values 1 + /. 
Calculated using Determinant QMC on 6 x 6 lattices with U — St and T — t/S. 
Error bars are within the symbol size when not shown. 



2.0 



1.0 



T = t/8 U = at 6x6 



f = 0.1 1 1 



-♦-f--f-#--f- 



Figure 5: Average density n on U = sites (triangles) and U = 8t sites (squares) 
separately for impurity concentration / = 0.111. The saturation value for n{U — 0) 
is about 1.71. Details as with Fig. 4. 
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T=t/8, U=8t, f=0.11 
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Figure 6: Average density n (filled symbols) and sublattice magnetization M (open 
symbols) on f/ = sites (triangles) and U = 8t sites (squares) separately. The 
saturation value for n{U = 0) is about 1.70. Calculated using DMFT for impurity 
concentration f = 0.11, U = 8t and T = t/8. Fig. 5 shows the corresponding results 
for n ior d — 2. 
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1 2 

Figure 7: Average total density n as function of fi for different impurity concen- 
trations / = 0, 0.056, 0.11, as calculated using DMFT. Dashed lines indicate the 
occurring values 1 + /. Details as with Fig. 6. Fig. 4 shows the corresponding results 
for d = 2. 
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one is not. This can be seen from the average sign which becomes very small, or from 
calculated quantities which occasionally deviate a lot from the majority of results 
(or even become unphysical, e.g. negative densities). To include such outliers in 
an average clearly makes no sense. A more reliable and robust way of determining 
the disorder average, i.e. one which is not very sensitive to the type of distribution 
which governs the outcomes, is the so-called midmean, or: 25% trimmed mean (an 
average over the middle half of the ordered data set |]2^). We have used this mean 
for the d = 2 QMC calculations; using another robust mean, the broadened median 
[P^, gives very similar results. We like to stress that this approach is not limited 
to the disordered case; in fact, we have applied it to the pure (/ = 0) model off 
half-filling as well (Fig. 4) by having different starting configurations take the role 
of different disorder realizations. 

In Fig. 4, it is possible to discern for non-zero / two flat regions in n(/i), i.e. 
two hardly compressible states. The first still occurs for n = 1, but has a much 
reduced charge gap, whereas a second gap develops at a density somewhat below 
n = 1 + f. One would expect a plateau exactly at n = 1 + / if the U = sites would 
fill independently of the U sites (because above 1 + / the U sites would have 
to start to be doubly occupied). Since the density at f/ = sites does not saturate 
at 2, but at a lower value (n = 1.76 and 1.71 for / = 0.056 and 0.111, respectively; 
see Fig. 5), clearly the two sets of sites are coupled. On the other hand, one notes 
from Fig. 4 that the second plateau for non-zero / terminates almost exactly at the 
l^c for / = 0. Therefore, the chemical potential needed to force double occupation 
of the U ^ sites is unaffected by the presence of f/ = impurities. The fact that 
n at the f/ 7^ sites remains pinned at 1 (Fig. 5) up to about the same value of fj, 
as for / = confirms this. As discussed above, indeed the error bars for non-zero / 
are smaller than for the pure case, indicating a less-severe sign problem. 

The remnant gap at n = 1 is in fact the more remarkable one, since the gap at 
= 1 + / is a clear descendant of the one at = 1 in the pure model (see above). 
This can also be appreciated from considering our model first with t = 0; for that 
case n{fi) will have two plateaus, at n = 1 -|- / and n = 1 — f (not mentioned above, 
but present because of particle- hole symmetry). Turning on the hopping not only 
leads to a rounding of the sharp steps, but also to the appearance of a new plateau 
at n = 1. We argue the gap at n = 1 to be due to induced AF order on the U = 
sites. To support this we have computed within DMFT for / = 0.11 besides the 
densities also the sublattice magnetizations M (Fig. 6). The decrease of M is clearly 
coupled to the pinning of n at 1, both for U = and U = 8t sites. For small n there 
is AF order on both types of sites, whereas for /i > 1.6 the density on U = sites 
saturates, the density on U = 8t sites starts to increase, and AF order disappears. 
Therefore, certainly for this rather large value of U, both charge gaps turn out to 
be intimately linked to AF magnetic order. 

Somewhat surprisingly, M on the U = sites becomes negative when n on 
these sites saturates. The reason is that electrons at f/ = sites are more strongly 
localized when their spin is parallel to the neighbors than when it is antiparallel 
(because of Pauli's principle). Therefore, the net moment on U = sites is parallel 
to its neighbors, i.e. opposite to the total staggered magnetization. 

For completeness we also show the total density as function of yU computed using 
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T=t/8 U = 8t 




n 



Figure 8: Kinetic energy | {K) \ as function of density n for impurity concentrations 
/ = (open triangles) and / = 0.11 (filled triangles). Calculated using Determinant 
QMC on 4 X 4 lattices with t/ = St and T = t/8. 



DMFT in the thermodynamic limit (Fig. 7). The overall picture of the d = 2 QMC 
calculations on 6 x 6 lattices (Fig. 4) is confirmed, although the plateaus close to 
1 + / are less well developed, certainly for / = 0.11. Another difference is that things 
happen at a slightly higher value of /x, because fewer fluctuations are accounted for 
in DMFT. 

Finally, we strengthen our case for an incompressible state close to n = 1 + / by 
showing in Fig. 8 the kinetic energy (as computed on a 4 x 4 lattice) as function of 
n. Clearly there is a dip in the electron mobility at n = 1 + /. A similar argument, 
using cusps in the kinetic energy as evidence for incompressible states, was found to 
work for the boson Hubbard model [^. We further note that additional support 
can be found in calculations of the density of states within DMFT that we present 



elsewhere 26 



In summary, using Determinant QMC calculations on finite lattices for d = 2 
combined with the infinite-dimension approach in the thermodynamic limit, we have 
studied the effect on ordered states in the Hubbard model of introducing disorder 
of the third kind, namely in the on-site interaction. We show that turning off the 
on-site interaction on a fraction of the sites of the lattice may enhance and stabilize 
AF order. It also leads to the occurrence of two incompressible states instead of one 
in the case without disorder; closer scrutiny of density and sublattice magnetization 
on the two types of sites teaches that both states are associated with AF order. 
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